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Double-diffusive instabilities are often invoked to explain enhanced transport in stably- 
stratified fluids. The most-studied natural manifestation of this process, fingering con- 
vection, commonly occurs in the ocean's thermocline and typically increases diapycnal 
mixing by two orders of magnitude over molecular diffusion. Fingering convection is also 
often associated with structures on much larger scales, such as thermohaline intrusions, 
gravity waves and thermohaline staircases. In this paper, we present an exhaustive study 
of the phenomenon from small to large scales. We perform the first three-dimensional 
simulations of the process at realistic values of the heat and salt diffusivities and provide 
accurate estimates of the induced turbulent transport. Our results are consistent with 
oceanic field measurements of diapycnal mixing in fingering regions. We then develop a 
generalized mean-field theory to study the stability of fingering systems to large-scale 
perturbations, using our calculated turbulent fluxes to parameterize small-scale trans- 
port. The theory recovers the intrusive instability, the collective instability, and the 7- 
instability as limiting cases. We find that the fastest-growing large-scale mode depends 
sensitively on the ratio of the background gradients of temperature and salinity (the 
density ratio). While only intrusive modes exist at high density ratios, the collective and 
7-instabilities dominate the system at the low density ratios where staircases are typically 
observed. We conclude by discussing our findings in the context of staircase formation 
theory. 
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1. Introduction 

When the density of a fluid depends on (at least) two components, nominally stably- 
stratified systems can, under certain circumstances, undergo double-diffusive instabilities 
leading to significant vertical buoyancy transport. Here, we focus on the case of the 
"fingering" instability, which often occurs in fluids which are thermally stably stratified, 
but have an inhomogeneous composition. A well-known example is found in upper layers 
of the Earth's oceans where evapo ration exceeds precipitation, lea ding to warm salty 
water overlaying colder fresh water (|Sternlll960l : ISchmitt et a/.ll2005l ). Since heat diffuses 



2 



A. Traxler, S. Stellmach, P. Garaud, T. Radko and N. Brummell 



faster than salt, parcels of fluid displaced downward rapidly lose their heat excess while 
maintaining a larger salt concentration. They become denser than the environment and 
continue to sink, forming structures called "salt-fingers" . Similar fingering instabilities 
can occur in any other thermally stably stratified solution, provided the con centration 
of the slower- diffusing solute increases with height ( Stern 1960; Sch mitt 1983 ). 

The saturated state of this instability, fingering convection, takes the form of ti ghtly- 
packed, vertically- elongated plumes of sinking dense fluid and rising light fluid ([Stern 
19601: iKunze l 1200 3 =). and significantly enhances the vertical transport of both heat and 



chemical composition. I n the ocean, fingering convection increases dia pycnal mixing 
within extended regions (ISchmittlll99l iKluikov fc Karlinlll995t IYou1I2002) of the thermo- 
cline by at least two orders of magnitude over molecular diff usion. It has been argued that 
the nutrient supply of the upper ocean (jDietze et al. the su rface temperature an d 

the surface fluxes of CO2 and O2 are all affected by this process (jGlessmer et aZ.1120081 ). 
Conditions fav orable for fingering convection also exist in many other natural systems 
( Schmittlll983l ). In the astrophysical context for example, a variety of situations lead to 
the development of unstable mean molecular weight gra dients in otherwise stably strat 



ified "radiative" regions within stars and giant planets ([Vauclairl l2004t IStancliffe et al 



2Q07t ICharbonnel fe Zahnl 120071 ). Since the long-term thermal evolution and chemical 



stratification of these objects is regulated by the transport bottleneck caused by radia- 
tive regions, the presence or absence of mixing by fingering convection can influence their 
observable properties dramatically. 

While fingering convection is by nature a small-scale phenomenon, it also has an in- 
triguing propensity to generate dynamical structures on very large scales, such as internal 
gravity wav es, th ermoha l ine in trusions and thermohaline staircases. As first argued by 
Stern (1969) and iHolve n (fl98l . a homogeneous field of fingers can become unstable to 
the so-called "collective instability" leading to the spontaneous generation of internal 
gravity waves in regions o f active salt fingering. This instability was later confirmed in 
laboratory experim ents by lStern fc Turner (1969) and in direct numerical simulations by 



Stern et a l. (2001). 



Thermohaline intrusions are different kinds of large-scale structures which are never- 
theless also often associated with fingering convection. They take the form of laterally 
interleaving layers with distinct temperature and salinity signatures, and can sponta- 
neously form in fluids which are s tratified both vertically an d horizontally. They are 
commonly observed in the ocean dRuddick fc Richards! l2003h. and have been repro- 
duced hi_l_ab_expj3™ 1999f ) and numeri- 
cally (|Simeonov fc Sternll2007l ). Seefauddick fc Kerrl 12003) for a detailed discussion of 
intrusion theory. 

Finally, one of the most dramatic signatures of active fingering convection in the ocean 
is the formation of mixed layers separated by salt finger interfaces, known as thermo- 
haline staircases. Persistent staircases have been documented in the Tyrrhen ian Sea, 
below the Mediterranean outflow, and in the western tropical No rth Atlantic (Schmitt 
1994. Layer format ion is also observed in laboratory experiments (jStern fc Turnerll 19691 



Krishnamurtill2003h. Layering enh ances vertical mixing by up to an order of magnitude 
(jSchmitt et aLll2005t IVeronial2007l ) relative to globally similarly stratified regions charac- 
terised by a smoother stratification. Fo rty years after the d iscove ry of this phenomenon 
in oceanographic field me asurem e nts (jTait fc Howd 1 19681 Il97ll ). a generally accepted 
explanation is still lacking. iRadkol (|2003l ) argued through theoretical and numerical anal- 
yses that the observed layering is likely to be caused by the so-called 7-instability — an 
instability driven by variations in the ratio of the turbulent heat and salt fluxes. 

The conventional approach to analysing the spontaneous generation of structures from 
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fingering convection uses the assumed separation of scale between the finger scale and the 
emerging structure scale to construct a mean- field theory, in which the effect of the small- 
scale fingering is modeled through turbulent fluxes. The resulting mean-field equations 
describe the evolution of the large-scale fields only and can straightforwardly be analysed 
for linear stability. Globally speaking, the modulation of the background stratification by 
large-scale temperature and salinity perturbations induces a modulation of the turbulent 
fluxes. When the divergence or convergence of these modified fluxes act to enhance the 
original perturbation, large-scale modes of instability are excited. Different variants of 
mean-field models have been individually successful in representing the gross properties 
of each of the aforementioned large-scale phenomena (intrusions, collective instability, 
7- instability). 

In this paper, we show that these various modes of instability can actually be de- 
scribed by a single unifying mean-field formalism, and are all recovered as limiting cases 
of our theory — each one corresponding to a different feedback mechanism between the 
large-scale perturbation and the induced turbulent fluxes. In Sj2] we present our unified 
mean-field model, and its relationship with previous work. In fJH we then perform a 
series of 3D simulations for parameter values typical of salty water in the ocean, de- 
signed to measure the turbulent transport of heat and salt as parametric functions of the 
background stratification. 

Using the small-scale flux laws derived, we then calculate and discuss in fj4]the expected 
growth rates of the various large-scale modes of instability as functions of the overall 
stratification of the region. Our results indicate that the relative importance of these 
various modes is highly sensitive to the density ratio (the ratio of the vertical temperature 
and salinity gradients normalised by their expansion/contraction coefficients). For low 
density ratio the dynamics of the system are primarily controlled by the collective and 
7-instabilities. For intermediate density ratios, the 7-instability is suppressed and the 
dynamics are dominated by gravity waves, with intrusive modes gaining importance. For 
larger values of the density ratio, only intrusive modes are unstable. 

Finally, we discuss our findings in SJSJ focussing on the implication of the measured 
turbulent fluxes for oceanic mixing in and on the role of the large-scale instabilities 
studied in the formation of thermohaline staircases in ^5.21 



2. Generalised mean-field theory of fingering convection 

2.1. The governing equations for homogeneous fingering convection 

Fingering convection, when observed in natural systems, typically occurs far from physi- 
cal boundaries. For this reason, we adopt an approach which minimises boundary effects 
by considering triply-periodic temperature, salinity and velocity perturbations driven by 
a steady and uniform fingering-unstable background stratifica tion. This setup has been 
advo cated by others before for studying fingering convection ( Stern et al. 120011 iRadkol 
2003), and is ideally suited to numerical simulations using spectral methods (see Sj3] and 



Paper II). It is important to note that it does not suffer from the well-known pathology 
of thermal conve ction in a triply-periodic system - the so-called homogeneous Rayleigh- 
Benard problem ( Borue fc Orszad[l996t ICalzavarini et aZ.ll2006l) . where the fastest grow- 
ing modes span the entire domain and depend sensitively on the aspect ratio of the box. 
Instead, the typical length scale of convective motions in the fingering regime is set by 
the diffusive length scales and is independent of the box size, provided the box is large 
enough (see Appendix). 

We consider a Cartesian coordinate system (x, y, z) with z increasing upward in the 
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vertical direction. In all that follows, we use the Boussinesq approximation. We assume 
that the background temperature and salinity profiles Tq{x,z) and So(x, z) are bilinear 
functions of x and z, Tq(x, z) = Tq x x + Tq z z and Sq(x, z) = Sq x x + Sq z z. Without loss 
of generality, the background fields are assumed to be two -dimensional (2D) by ali gning 
the horizontal gradients with the x-axis. We assume (as in lWalsh fc Ruddicklll995h that 



the overall horizontal density gradient is zero, in which case aT 0x = (3Sq x where a and j3 
are the coefficients of thermal expansion and compositional contraction respectively. The 
slope of the background temperature gradient in this coordinate system is <fi = Tq x /Tq z . 

We perform a standard non-dimensionalisation pro cedure for s tudying local finger- 
ing convection. We use the expected finger scale (see Stern 1960l ) as the length scale, 



d = (k,tv I gaToz) 1 / 4 , where g is gravity, v is viscosity and kt is thermal dif- 
fusivity. We then define the corresponding thermal diffusion time scale, [t] — d 2 /Kr, 
the velocity scale [u] — Kx/d, and the temperature and salinity scales, [T] = T$ z d 
and [S] — (a/ f3)To z d. Nondimcnsional parameters of interest are the Prandtl number, 
Pr = v/kt, the background density ratio, Rq = qTq z / (3Sq z and the diffusivity ratio, 

T = K S /kT- 

The non-dimensional equations for the evolution of the velocity field u = (u, v, w) and 
the temperature and salinity perturbations T(x,y,z,t) and S(x,y,z,t) arc then: 

'" I — +u- Vu] = -Vp+(T-S)fc + V a u, (2.1a) 



Pr V dt 

V-u = 0, (2.16) 

dT 

— +<j)u + w + u-\7T = V 2 T, (2.1c) 

f)Q 1 

— +<f>u+ — w + u-VS = tV 2 S, (2Ad) 

Ot Rq 

where p is the non-dimensional pressure perturbation from hydrostatic equilibrium and 
k is the unit vector in the z-direction. 

2.2. Generalised mean-field theory 

As discussed in SJTJ fingering convection is often associated with the emergence of dy- 
namical structures on scales much larger than individual fingers. We begin by deriving 
a generalised set of mean-field equations, and then study their linear stability to various 
large-scale modes. 

2.2.1. Mean-field equations 

As in iRadkol t|2003h . we are interested in the large-scale behaviour of the system of 



equations (I2.1aj|2.1d| when averaged over spatial/temporal scales of many fingers. We 
introduce the notation 7TT , where the overbar denotes an averaging process which we 
assume may commute with spatial and temporal derivatives. Let u = u+u' and similarly 
for T and S, in which case u' = T' = S' = 0. The averaged governing equations now 
become 

-j- f^+u-Vu) =-Vp+(T-£)k + V 2 u--|-V-R, (2.2a) 
Pr \ ot I Pr 



dT 
~dt 



VT = V T — V • Ft, (2.26) 



dS 1 — — 

— + (jm+ — w + u- VS = tV 2 S - V ■ F s , (2.2c) 

Ot Rq 
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where Rij — u^u^, and the turbulent fluxes are Ft = u'T', F$ = u'S'. 

In what follows we now drop the overbar and only refer to the evolution of the large- 
scale fields u, T, S. As in previous analyses, we assume that the Reynolds stress term 
is small enough to neglect (a fact that is easily verified a posteriori), and additionally 
assume that only the vertical component of the heat and salt fluxes are large enough to 
be significant (so that Ft ~ -Prk, F$ ss -Flsk). However, we retain the diffusion terms in 
all three equations. In nondimensional terms, the turbulent fluxes are characterised by 
the Nusselt number Nu and the turbulent flux ratio 7, defined as: 

-(i + dT/dz) ' {2 - 6) 

7=f- (2-4) 



Note that these definitions of Nu and 7 differ somewhat from those of Radko (2003), who 



includes the molecular diffusive terms in his definition of Ft and Fg- Our formalism has 
greater generality, since it allows for horizontal diffusive fluxes. The difference becomes 
important at low Nusselt number. 

We now make the key assumption that at any given time both Nu and 7 depend only 
on the local value of the density ratio R p , which in nondimensional terms is 



aT 0z (l + dT/dz) 
"~ pS 0z {l + (%fc)BS/dz] 



Rn — 



I + dT/dz 

- R °l + R a dS/dz- (2 - 5) 

The functions Nu(i? p ) and j(R p ) can be determined experimentally, using numerical 
simulations (see §5$. 

The system of equations describing the evolution of the large-scale quantities u, T, S 
is now 

^ f^+u-Vu) =-Vp+(T-,S)k + V 2 u, (2.6a) 

— +4>u + w + u-VT = V 2 T- — -, (2.66) 
at oz 

^- + ( / ) u + ^w + u-\7S = tV 2 S~^, (2.6c) 
at Rq oz 

where the flux derivative terms on the right hand side use (|2.3p and (|2.4p to express Ft 
and Fg in terms of Nu and 7. 



2.2.2. Linearised mean-field theory 

The mean-field equations derived above exhibit steady solutions describing a state of 
homogeneous fingering convection, with zero mean velocity, zero deviation from the back- 
ground temperature and salinity fields and constant (non-dimensional) heat and salinity 
fluxes Fto = (1 — Nuo(-Ro)) and F$o — j(Ro)/Fto- The stability of this homogenous tur- 
bulent state can be investigated by adding a small perturbation to the mean quantities, 
and linearising the mean-field equations. Large-scale temperature and salinity perturba- 
tions induce large-scale variations in the density ratio, so that R p = Rq + R' p . This, in 
turn, modulates the turbulent fluxes in a way which may in some circumstances further 
enhance the initial perturbations or quench them. Various modes of instability are re- 
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lated to different feedback mechanisms between the fields and turbulent fluxes, through 
the parametric functions Nu(i? p ) and j(R p ). 

Assuming that perturbations away from the linear background gradients are small, we 
first have, expanding (|2.5[) to linear order: 



(2.7) 



R p — Ro + -g^ — Ro~q~*J — Ro + R' p - 



which uniquely defines R' p , and then 



Nu(i? p ) « Nu(flo) 



dNu 



dR n 



R' 



(2.8) 



Ko 



and similarly for 7. Rearranging (|2.3|) and (|2.4I) yields Fx = (1 — Nu)(l + dT/dz) and 
F5 = Ft/ j- It then follows that 



dz 
dF s 
dz 



8 2 T 



d 2 S 



^1(^2 -^o^-j )(Nuo-l) 



d 2 T 
dz 2 
1 <9F T 

7o 9z 



9z 2 <9z z 

where we have abbreviated Nu(-Ro) = Nuo, j(Ro) = 70, an d defined 

97- 1 



Ai = R 



A 2 = R 



dR P ' 

dm 



dR n 



(2.9) 
(2.10) 

(2.H) 
(2-12) 



Note that our A\, A 2 are not strictly equal to those defined by iRadkol ( 2003 ) but reduce 
to the same quantities in the limit where turbulent fluxes are much larger than diffusive 
fluxes. 

A standard linear stability analysis of (|2.6p . using normal modes of the form {u, T, S} — 
{u, T, 8} exp (At + ilx + imy + ikz), yields a cubic equation for the growth rate, A 3 + 
a 2 A 2 + aiA + ao = 0, with 



a 2 = |k| 2 (I + Pr + r) + fc 2 
ai = |k| 4 (rPr + r + Pr) + fc 2 |k| 2 

- J 4i J R (l + Pr)(Nu -l 
a = |k| 6 rPr + fc 2 |k| 4 Pr 



(I-A 1 i? )(Nuo-l) + A 2 1- 



Rp 

Jo 



(2.13a) 



(r + Pr)(A 2 + Nu - 1) - A 2 (l + Pr)^ 

7o 



- £r%i? (Nu - l) 2 +Pr 



I 2 l 2 + m 2 



1-^), (2-136) 



(r-A 1 i? )(Nuo-l) + A 2 (r-^ 

7o 



+Pr- 



I 



I 2 + m 2 



+fc 2 A 1 (l- J R )(Nu -l)(/^ k<j>) 



I 2 



-k 2 [A 2 (l-i? )+Nu -l] 



I 2 + m 2 



1 



1 



-Ro 7o 



k(j> 1 



fc 4 |k| 2 Pri?o-4i(Nu - l) 2 



(2.13c) 



To 



where Ikl 2 = k 2 + I 2 + m 2 . 
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2.2.3. Relationship with previous theories 

Presented above is a unified formulation of several mean-field theories, including the 
effects of all diffusion terms and the contribution of variable turbulent flux ratio 7 as well 
as allowing for the presence of lateral background gradients in temperature and salinity. 
Several limiting cases have been discussed previously in the literature (unless otherwise 
noted these are 2D theories, so m = 0). 

The fingering instability. Although technically not a mean- field instab ility, it is re- 
assuring to note that the fingering instability itself (e.g. iBaines fc Gilllll969h is recovered 
when turbulent fluxes and lateral gradients are ignored (A\ = A2 = Nuo — 1 = </> = 0). 
In that case the cubic defined by (|2.13p becomes the well-known cubic equation for the 
growth rates of the fingering modes with 

a 2 = |k| 2 (l + Pr + r), (2.14a) 
a x = |k| 4 (rPr + r + Pr) + Pr^ (l - i-) , (2.146) 

a = |k| 6 rPr + PrZ 2 ( r - — ) . (2.14c) 



The collective instability, as derived bv lStern et al. ( 2001 ). is recovered from (|2.13[) 
by omitting lateral gradients (0 = 0), neglecting possible variation in 7 (so A± = 0) and 
discarding the diffusion terms for temperature and salinity — but not velocity — in (|2.2|) . 



a 2 = Pr|k| 2 + fc 2 



ai = Prfc 2 |k| 2 A 2 



An 



^ Ro 

7o 



Nu - 1 



- Ro 

7o 



Nu - 1 



Pr- 



1 

Ro 



a 



k 2 ! 2 



R ) + Nu ~ 1] 



(2.15a) 
(2.156) 
(2.15c) 



1 1 

70 -Ro, 

Stern ( 19691 ) argued that modes are excited when the Stern number, written in our 



notation as 



A = 



(Nu - 1) 



1 



Pr 



1 Ro 



(2.16) 



exceeds a value of order one. The unstable modes essentially represent overstable gravity 

waves. _ ^ 

\ 1 elegant physical interpretation ( Stern et al. 200lh of the collective instability is 
obtained by analogy with the laminar, linear double-diffusive instability in the "diffusive 
regime", where the slowly diffusing field is stably stratified while the rapidly diffusing 
field is unstably stratified. In this case, growing oscillatory modes akin to internal gravity 
waves are excited instead of fingers. Since fingering convection induces a mean salt flux 
larger than the heat flux the roles of the two fields are reversed, and the faster diffusing 
field is now the salinity field. From a turbulent point of view, the "diffusive regime" is 
recovered. 



The theory of intrusions of Walsh fc Ruddickl ( 19951 ) is recovered by discarding the 
diffusion terms in (|2.2j) and setting 7 constant (A\ =0), as well as neglecting Reynolds 
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stresses in their formulation: 



<>> = ir A 2 (Nun - 1)(1 - AiRo) 

I 2 ( 1 

a x = -fc 4 i?o(Nu - l) 2 Ai + Pr— 1 - — 



Ro 



Ik 2 ( 

+A 1 (l - i? )(Nu - - H)} . 



1 

7o 



' ' 7o RO 



(2.17a) 
(2.176) 
(2.17c) 



The mechanism underlying intrusive instabilities can be illustrated by imagining an al- 
ternating horizontal s hear flow superimposed on the background stratification (see e.g. 
Ruddick fc KeTrll2003h . The lack of horizontal density variation in the background (where 
lateral gradients of temperature and salinity compensate) implies that the background 
isohalines are steeper than the isothermals, and are therefore more strongly affected by 
horizontal advection. Alternating vertical variations in R p result, which in turn strength- 
ens or weakens the fingering action, and the resulting flux convergences and divergences 
reinforce the intrusive motion. Dependin g on the orientation of t he perturbation, both 
direct and oscillatory modes are possible (jWalsh fc Ruddicklll995l ). 



The 7-instability, as derived in Radko ( 2003lh is recovered by considering horizon- 
tally invariant perturbations (I = m = 0, |k| 2 = k 2 ) with zero velocity field. From the 
remaining temperature and salinity equations we obtain a quadratic^ expression for the 
growth rate: 



0,2 = 1, 
ai = k 2 

a = k 4 



Rg 

7o 



1 + r + (1 - Ai#o)(Nu - 1) + A 2 I I 
(r - A 1 i? )(Nu - 1) - AMNuo - if + A 2 



Ro 

7o 



(2.18a) 
(2.186) 

(2.18c) 



Differences between these coefficients and those given in iRadkd (l2003h arise from our 
alternate definition of 7, but can be shown to be reduce to each other in the limit of 
large Nusselt number. Note however that for 7-modes, which do not have any horizontal 
variation, the use of the total fluxes F*, ot = F T - (1 + dT/dz) and i^ ot = F s - t(1/R + 
dS/dz) (originally advocated by Radko) recovers his much simpler quadratic with 



«2 

a 1 



I I ^'"./.'ulXUn-i- .1, (1-^ 

7o 



a = -^A^RoNul, 



(2.19a) 
(2.196) 
(2.19c) 



where 7^ = F*, 7F^ ot and A\ ot = i? d(l/7 tot )/dR p . As shown by Radko, a sufficient 
condition for the existence of a positive real root is that A^ > 0, or in other words that 
7 tot should be a decreasing function of R p . The physical interpretation of this s o-called 
"7-in stability" is fairly subtle, and is described in detail in the original paper (jRadkol 
200.^ . 



t To see why our formalism yields a cubic while Radko's yields a quadratic, it should be 
noted that (|2.13[) can be factored (with I = m = 0) as (A + fc 2 Pr)(A 2 + 61A + bo). The first root 
describes the viscous decay of any initial (horizontally invariant) velocity perturbation. 
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3. Turbulent flux laws 

To proceed forward and estimate growth time scales for the various mean-field modes 
of instability excited by fingering convection, we need to determine the non-dimensional 
turbulent fluxes Nu and 7 as functions of the density ratio. Naturally, these depend on 
the diffusivity ratios Pr and r relevant of the system studied. Here, we choose to focus 
on the case of salty water (Pr = 7 and r = 0.01), as it is directly applicable to the 
oceanographic context. Appendix A summarises the numerical algorithm, describes the 
experimental protocol for determining heat and salt fluxes, and discusses the problem 
of selection of the domain size to be used for these experiments. What follows are the 
results of a body of simulations at different density ratios R p . 

3.1. Typical results 

Figure [T] shows a visualisation of the salinity field obtained in the saturated state of a 
fingering system with R p — 1.2, R p = 2 and R p = 10. On account of the small diffusivity 
of salt compared with all other fields, a broad range of spatial scales exists, and a high 
numerical resolution is a priori required to resolve all arising structures and correctly 
model the system. We find that, as expected, the salinity field has a complicated struc- 
ture for small R p , but successively becomes more organised with increasing density ratio. 
Regular, vertically elongated filamentary structures dominate for R p 10. In order to 
ensure that the smallest scales of the salinity field are fully resolved, we had to use the 
highest resolution available (a grid of 768 x 768 x 1536) at R p — 1.2, although half that 
resolution is sufficient for R p ^ 2. Furthermore, it turns out that a rough estimate of the 
flux laws can in fact be made with a much coarser resolution (about 32 3 ). This rather 
surprising result shows that the diffusion of salt does not play an important role in con- 
trolling the mixing in the heat-salt system, for very turbulent flows (low R p ). Moreover, 
it suggests that low-resolution simulations may be sufficient to estimate turbulent fluxes 
for any high-Pr, low-r fluid. 

3.2. Turbulent flux laws for the heat- salt system 

The control parameters used in each simulation, along with some key results, are sum- 
marised in table [TJ Plots displaying the most important findings are shown in figure [H 
which also contains results from an accompanying set of 2D simulations. 

As expected, we find that the turbulent fluxes \Fs\ and \Ft\ decrease rapidly with 
increasing density ratio, and tend to be considerably larger in the three-dimensional 
(3D) case than in 2D. The ratio of 3D to 2D fluxes is not constant, but tends to grow 
with increasing R p . Because of their obvious oceanic relevance, tabled] also lists values for 
the "eddy" diffusivities Kt — kt\Ft\ and K$ — ktR p \Fs\ of heat and salt. At R p = 1.2, 
we find Kt ~ 0.21 cm 2 /s and K$ ~ 41 cm 2 /s. Both quantities quickly decrease with 
increasing R p . 

The turbulent flux ratio 7 tends to be larger in 2D than in 3D. It initially decreases 
quickly with growing R p , attains a minimum around R p « 7 (7 « 0.45 in 3D) and then 
slowly increases again. A widely used theoretical pr ediction based on linearly fastest 
growing modes, originally proposed by ISchmittl (|l979t ). tends to over-estimate 7 consid- 
erably, and also predicts the minimum to occur at a smaller (R p = 4) value of R p . The 
total flux ratio 7 tot (see §2.2|) , which plays a prominent role in the 7-instability theory, is 
shown in figure El Its value deviates significantly from the turbulent flux ratio 7 at higher 
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Figure 1 . Snapshots of the salinity field S in simulations of fingering convection in the heat-salt 
system (Pr — 7, r = 0.01). a) Salinity field at R p — 1.2, plotted on the three planes x = 0, y — 
and z — L z . b-d) Volume rendering of the salinity field for R p = 1.2, R p = 2 and R p = 10 (from 
left to right). In all cases, the simulation domain contains 5 x 5 x 10 FGW (see main text). Note 
how the typical amplitude of the salinity perturbation in a finger is of the order of 1/R p t, or, 
in dimensional terms, dSoz/f. 



values of R p , where the Nusselt number is lower. As a result, the position of the minimum 
of the curve occurs for lower R p , thus restricting the range for which the 7- instability is 
expected to R p ^ 4. Furthermore, since the growth rate of the instability is proportional 
to d(l/7 tot )/dR p , we find that 7-modes should only be significant for R p < 2. 

Finally, we find that the Stern number A, which controls the dynamics of the collective 
instability, exceeds unity for R p 7 in the 3D case, while two-dimensional simulations 
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Figure 2. Parametric dependence of the non-dimensional fluxes Ft, Fs as well as their ratio 7 
and of the Stern number A as a function of R p . Results from both three- and two-dimensional 
simulations are shown. Panel c ) also cont a ins a theoretical prediction of j(Rp) based on the 
fastest growing linear modes (Sch mittlll979l '). 





R p = 1.2 


Rp = 1.5 


R p = 2.0 


Rp =4 


R P = 7 


Rp = 10 


resolution 


768 2 x 1536 


768 2 x 1536 


384 2 x 768 


384 2 x 768 


384 2 x 768 


384 2 x 


768 


Afaverage 


39.1 


57.8 


121.2 


223.8 


422.7 


390.1 




\F T \ 


153.5 ±11.7 


73.2 ±5.7 


37.6 ± 2.2 


13.3 ± 1.0 


5.48 ± 0.38 


3.35 ± 


0.21 


m 


241.8 ± 13.1 


126.4 ± 7.5 


70.3 ±3.1 


27.4 ± 1.5 


12.1 ± 0.65 


7.29 ± 


0.34 




0.63 ±0.02 


0.58 ±0.01 


0.53 ±0.01 


0.49 ±0.01 


0.45 ±0.01 


0.46 ± 


0.01 




14.1 ±0.4 


9.4 ±0.3 


6.5 ±0.15 


3.7 ±0.1 


2.37 ±0.06 


1.82 ± 


0.04 


K T [10- 6 m 2 /s] 


21 ±2 


10 ±1 


5.3 ±0.3 


1.9 ±0.1 


0.77 ±0.05 


0.47 ± 


0.03 


K s [10- 6 m 2 /s] 


41 ±2 


27 ±2 


20 ±1 


15 ± 1 


11 ±1 


10 ±0 


5 


A 


76 ±3 


23 ±1 


9.4 ±0.3 


2.7 ±0.1 


1.1 ±0.05 


0.63 ± 


0.03 



Table 1. Summary of simulations of fingering convection for the heat-salt system 
(r = 0.01, Pr = 7) in a computation al domain containing 5 x 5 x 10 fastest growing finger 
wavelengths (FGW, see I Schmittl 1 1 9 791 ) . Afaverage denotes the length of the time interval over 
which the data has been averaged, Kt = ktI-FtJ an d Ks = KtR p \Fs\ are the "eddy" diffu- 
sivities for heat and salt, taking kt = 1.4 x 10~ 7 m 2 /s, and A is the Stern number defined in 
p36jl . 
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■ Total 

O Turbulent only 



Figure 3. Comparison of the flux ratio 7 tot as calculated using the total heat and salt fluxes 
(both turbulent and convective), and 7 as calculated using only the fluxes due to turbulent 
fingering convection. The diffusive contributions become important as R p increases and the 
fingering fluxes drop, affecting not only the values of the flux ratio but also the location of the 
minimum of the curve. 



considerably under-estimate A and therefore underestimate the range of R p for which 
the system may be unstable to these modes. 

Further discussion of the implications of these simulations is deferred to For now, 
the above results provide the necessary data to apply the mean- field theory of |}5]to the 
oceanic parameter regime. 



4. Dominant modes of instability as a function of background density 
ratio 

The flux laws determined above enable us to estimate the growth rates for the various 
mean-field modes of instability discussed in [21 As seen in i )2.2.31 up to four modes of 
instability exist, but it is not immediately clear which mode dominates in the various 
regions of parameter space. To answer this question in the oceanic context, we now 
examine the solutions of the growth rate equation (|2.13[) for different values of Rq , with 
the corresponding 70, Nuo, A% and A2 calculated from the turbulent fluxes measured 
in section |3] and shown in table 2. We then find the largest growth rate for a given 
mode geometry (as determined by k and I), maximising Re(A) over the three roots of 
the cubic. Figure [4] shows this maximum growth rate as a function of wavenumber for 
three representative values of the background density ratio: i?o = 7, where the fingering 
instability is weak, R = 1.5, where the density gradient is close to unstable and turbulent 
fluxes are large (see tabled, and an intermediate value of Rq = 4. The plots show I on a 
logarithmic scale to capture the wide range of relevant horizontal lengths, from the small 
filaments of the fingering instability up through extensive lateral intrusions. 

Since the growth rate of the fingering instability is recovered from our mean-field theory 
when Ax, A 2 , etc. are zero (as discussed in H'2.2. 3[) . an analogous feature appears here 
even though Ai,A 2 ^ 0. This "fingering" mode appears as a "bulb" on all plots, at 
the smallest horizontal scales (large I) and large vertical scales (low k). Note, however, 
that mean-field theory should not be applied to model such small-scale structures. In 
practice, the bulb merely serves to indicate the region of I space (log/ —1) above 
which the theory is no longer applicable. 
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Table 2. Empirically derived parameters for growth rate prediction using the unified theory at 
high, midrange, and low values of background density i?o- Ai and Ai are estimated by using 
neighboring values of the density ratio. 

R Nu 70 A 1 A 2 



7.0 5.48 0.45 0.16 -11.61 
4.0 13.3 0.49 0.27 -25.7 
1.5 73.2 0.58 0.56 -217.31 



At high density ratio (top of figure [4]), in the absence of lateral gradients, only the 
fingering mode remains. The presence of a lateral gradient introduces two additional 
regions of instability, one oscillatory and one direct, both confined to large vertical and 
horizontal scales (i.e., small k, I). As expected, lateral gradients break the symmetry of 
the solutions, since the k(f> terms in (|2.13p distinguish between positive and negative k 
perturbations. The direct mode corresponds to a lateral intrusion which typically grows 
on a time scale of about 30 hours, with a horizontal scale of the order of a kilometer and 
vertical scale of a few metres (i.e., with a slope of the order of (f). 

For an intermediate density ratio (Rq = 4, middle panels of the figure), the gravity 
waves of the collective instability appear at a range of vertical and horizontal scales 
starting at I = 0.055, k = 0.06 (a physical size of about a metre in each direction), with 
a growth time scale of about 30 hours. The lateral gradient strongly modifies the gravity 
waves, increasing both their maximum growth rate and the size of the instability region 
for negative k values, while suppressing growth for positive k. As with Rq = 7, the lateral 
gradient also triggers a direct intrusive mode at large horizontal scales. 

Finally, at low density ratios the system is dominated by the collective instability 
(oscillatory) and the 7- (direct) instability, and is now unstable to a continuous range of 
modes on both large and small horizontal and vertical scales. At the scales for which the 
mean- field theory is valid (k <C 1), the collective instability grows fastest, most unstable 
on scales of a metre both horizontally and vertically, and with a growth time scale around 
two hours regardless of the presence of a lateral gradient (0 = 0.01). For comparison, a 
7-mode of the same vertical scale grows at roughly a third this rate. 

Figure [5] shows the growth rate of various modes in our theory at low density ratio 
(R p = 1.5) and in the presence of lateral gradients {4> = 0.0 1). Also shown are the 
growth rates of the collective instability of IStern et all <l200lh . of th e 7-instability of 
iRadkol (2003), and of the intrusive instability of Walsh fc Ruddickl ( 1995 . absent Reynolds 
stresses). For a horizontally invariant perturbation (I — 0, equivalent to a domain width 
of 2ir/l — s- 00), the mean-field theory matches Radko's 7-instability. For perturbations 
with moderate horizontal scales (2ir/l = 200d, approximately two metres in dimensional 
terms), the unified theory predicts the presence of gravity waves with similar vertical 
scales and recovers the collective instability. As the vertical wavelength decreases, the 
mode gets flatter, and becomes a 7-mode. At the largest horizontal scales (2n/l = 2-10 5 d, 
about two kilome t ers), our theory recovers the growth rate of intrusive modes from 



Walsh fc Ruddickl (119951) . 



Using this unified formalism, we have therefore demonstrated how the dominant type of 
large-scale instability in salt fingering systems depends strongly, not only on the scale of 
the perturb ations considered, but also on the background density ratio. Field observations 
l Youll2002h reveal that R p can vary significantly in the ocean. In nearly unstable regions 
(R p — > 1), the collective and 7-instabilities control the dynamics of the system (also see 
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-0.05 0.05 -0.05 0.05 

k k 

Figure 4. Predicted real part of growth rates for the fastest growing perturbations, where 
colour is scaled to logRe(A). Only positive values are shown; the region is left white if no modes 
grow. The horizontal axis shows vertical wavenumber, and the vertical axis shows the logarithm 
of horizontal wavenumber to capture the broad range of expected scales. This particular dis- 
play choice yields these characteristic "flower-plots." The left-hand column shows results in the 
absence of lateral gradients (cj> = 0) while the right-hand column shows results for a typical 
oceanic value of </> = 0.01. In each of the six figures, regions surrounded by a dark contour show 
oscillatory behaviour, by contrast with the direct modes. For example, the symmetric "bulbs" 
at high / are direct modes, corresponding to the growth rate of individual fingers. Top: High 
density ratio (_Ro = 7). Middle: Midrange density ratio (Ro = 4). Bottom: Low density ratio 
(Bo = 1.5). 
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Figure 5. Comparison of the real component of projected growth rates between theories at 
Pr = 7, r = 1/100, for R = 1.5 and <j> = 0.01. The I = ( 2n/l = oo) mode growth rate is 
identical to the one predicted by the 7- instability of Radko (2003); the remaining lines have 
horizontal perturbation wavelengths of 200d = 1.8 m and 200000c! = 1.8 km. The mean-field 
theory captures both the collective instability at large vertical scales and the 7-instability at 
small vertical scales, and matches well with the intrusive instability at the largest scales. 

Paper II), but as the system becomes more strongly stratified, these modes disappear and 
lateral gradients emerge as the dominant factor in the creation of large-scale structures. 



5. Discussion and Conclusions 

5.1. Turbulent flux laws 

The simulations of H3.2I are of direct relevance to the problem of parameterising doublc- 
diff usive mixing in the ocean. These problems ar i se in basin-scale ocean circulation mod- 
els (|Gargett fc Hollowavlll993 IZhang et aZj|l998HMerrvneld et a?Jll999h and in fine-scale 
studie s focusing on the dynamics of intrusions, interna l waves, and thermohaline stair- 
cases (| Walsh fc Ruddickll200d : ISimeonov fc SternlEioi IStern fe Simeonovll20ol iRadkol 
2005). While several attempts ha ve already been made to deduce the small domain flux 



laws from n umerical simulations (|Shenlll995[ IStern et aill2001t IStern fc Simeonovll2005t 
Radkoll2008l) . the computational restrictions in early studies precluded direct treatment of 
the oceanographic case (characterised by three-dimensional dynamics, Pr=7, r = 0.01). 
Instead, simulations were either two-dimensional or employed diffusivity ratios signifi- 
cantly higher than the heat-salt value of 0.01. The double-diffusive flux laws were de- 
duced by extrapolation of the numerical results obtained in the computationally accessi- 
ble regime — an approach clearly requiring a posteriori validation. The simulations sum- 
marised in table[T]are the first DNS that meet the challenge of solving the actual heat/salt 
problem in three dimensions. A comparison of these results with earlier studies reveals a 
good qualitative agreement with earlier estimates. In re trospect, it is perhaps surprising 
to see how well the former educated guesses of flux laws ( Schmitt ll979HStern e7olll200lh 
captured the pattern of heat/salt diffusivities as a function of density ratio. 

The comparison with oceanographic field measurements is more ambiguous since small- 
scale mixing in the ocean is driven by a combination of double-diffusion and turbulence — 

their relative contribution is uncertain and much debated. Nevertheless, t he careful analy- 

sis of the NATRE (North Atlantic Tracer Release Experiment) data set bv lSt. Laurent fc Schmitt] 
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(1999) made it possible to evaluate the salt finger diffusivities directly from observations. 
The estimated salt diffusivity is characterised by a monotonically decreasing dependence 
on densit y ratio, reducing from Kg = 5 0-10~ 6 m 2 /s at R p = 1.4 to Ks = 10-10 -6 m 2 /s at 



R p = 1.8 ( St. Laurent &: Schmitt 19991 ). Once again, these values are in broad agreement 



with the DNS results summarised in table 1. 

Finally, the presented synthetic data (table 1) make it possible to assess the relev ance of 
severa l hypotheses proposed to explain the physics of equilibration. Most notably, ISternl 
(1969) suggested that the amplitude of salt fingers could be limited by the collective insta- 
bility, with equ ilibrium fluxe s characterised by values of A ~ 1. A similar suggestion was 
put forward byE unze (Il987l) . who pointed out that Stern 's (Il969l) criterion is equivalent 
to specifying the Richardson number based on scales of individual fingers and speculated 
that an increase in A above unity would be followed by the ra pid destructi o n of f ingers 
by secondary instabilities — an idea most recently revisited by llnoue et al. (2008). Our 



results, which reveal variation in A by two orders of magnitude, emphasise the limitations 
of Stern/Kunze hypothesis and motivate the search for alternative conceptual models. 

5.2. Large-scale instabilities and implications for the formation of thermohaline 

staircases 

In this paper, we have treated three proposed mechanisms for large-scale instability in 
salt fingerin g systems, unifying under one framework what had pre viously been studied 
in isolation (jWalsh fc Ruddicklll995t IStern et a?]|200lURadkoll2003l ). Note that a related 
approach to the problem, considering the effect of the Richardson number on the fluxes 
and focussing on thermohaline interleaving rather th an the oscillatory modes of the 
collective instability, may be found in the recent work of lSmvth fc Ruddicld (|201(lh . In our 
work, as in theirs, considering all instability mechanisms in a single formalism opens the 
possibility of comparing the growth rates of the various mean-field modes to one another 
and establishing the dominant ones in each region of parameter space. Furthermore, 
using the realistic flux laws discussed above, we are now able to give robust quantitative 
predictions for the presence and growth rates of each mode individually. 

For the heat-salt system, we find that no single one of the proposed instability mech- 
anisms is expected to dominate in all fingering-unstable regions of the ocean. At high 
density ratio (R p ^ 7), for example, the Stern number drops below one and only lateral 
intrusions may be destabilised. As shown in |j4j intrusive modes with horizontal scales 
on the order of a kilometer and vertical scales of a few metres are expected to grow on 
a time scale of about 30 hours. As the density ratio decreases below seven, gravity-wave 
modes are also destabilised, growing fastest at horizontal and vertical scales of a few 
metres. The relative growth rates of the gravity waves and the intrusive modes depend 
sensitively on their spatial extent and on the slope of the isothermal contours (</>) with 
respect to the vertical, in a manner which can be evaluated through our theory. Finally, 
when the background stratification is close to neutral stability, which is the case for most 
fingering regions of the ocean, 7-modes are also unstable and grow on similar time scales 
as the gravity waves — on the order of a few hours, much faster than the intrusive modes. 

These findings also enable us to place constraints on existing theories for the forma- 
tion of thermohaline staircases. Indeed, all three mean-field modes of instability have 
been proposed to generate these structures in the proce ss of their nonlinear development 
(|Walsh fc Ruddicklll995t IStern et aZJ[200lt lRadkoll2003h . However, it is important to note 



that staircas es are typically only observed to exist in very low density ratio environments 
(R p < 2, see ISchmittl 1 1 98 II ) . We find that in this parameter regime (see figure [5]), intru- 



sions grow one to two orders of magnitude slower than gravity waves or 7-modes unless 
lateral gradients are exceptionally strong (which could happen in some regions of the 
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Mediterranean outflow, for example). This limits the relevance of intrusive modes when 
applied to the formation of staircases in the bulk of the thermocline. We also find that 
gravity-wave modes grow faster than 7-modes on all spatial scales for which mean-field 
theory is applicable (vertical wavelength greater than about lOOd, see Appendix). This 
should a priori point to the collective instability as the mechanism responsible for layer 
formation. 

However, the aforementioned correspondence between the locations of observed oceanic 
staircases (R p < 2) and interval of strongly decreasing 7 tot (i? (:) ) is too remarkable to be 
dismissed. In addition, the onl y existing num erical simulation to date for which staircase 
formation has been observed ( Radko 20031 ) has unambiguously identified a 7- mode as 



the staircase precursor. For these reasons, the 7-instability could prove to be just as 
important as a generating mechanism for these large-scale structures. At this point, large- 
scale numerical simulations are the only avenue towards further progress, an avenue we 
follow in part II of this paper. 
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Appendix A. Numerical determination of flux laws 

A.l. Description of the numerical algorithm 

We measure the nondimensional turbulent fluxes of heat and salt as functions of the 
density ratio using the following numerical algorithm. We solve the original set of equa- 
tions (|2.1al [2~Td| for homogeneous fingering convection using, as explained in §2.11 triply- 
periodic boundary conditions for all perturbations, e.g. 

T(x, y, z, t) = T(x + L x , y, z, t) = T(x, y + L y , z, t) = T(x, y,z + L z , t) , (A 1) 

where {L x , L y , L z ) defines the dimensions of the computational box (in units of d). Note 
that in these units the global Rayleigh number of a simulation is equal to L\. In this 
section, all quantities refer to the full field containing all scales (by contrast with £j2T2j - 



We use a spectral algorithm based on the classical Patterson-Orzag method (jCanuto et al. 



20071) widely used for simulations of homogenous turbulence. Nonlinear products are eval- 



uated on a grid in physical space, and the 3/2-rule is used to avoid aliasing errors (for 
reference, N grid points in a coordinate direction corresponds to Fourier modes up to 
wavenumbers of (2/3)iV in that direction). Note that our simulation is a Direct Numeri- 
cal Simulation (DNS), with no subgrid scale model. A third order, semi-implicit Adams- 
Bashforth / Backward-Differencing algorithm ( Peyretll2002 ) is used for time-stepping. All 



diffusive terms are treated implicitly, while the advection terms are explicitly treated. The 
above time-stepping method was chosen since it offers a relatively large stability domain 
that includes a part of the imaginary axis at a comparatively low cost. In order to guar- 
antee numerical stability, the time step is adjusted dynamically. The code was designed 
to run efficiently on massivel y-parallel supercomputers and employs a transpose-based 
parallel transform algorithm (jStellmach fc Hansen! 120081 ). 
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A. 2. Experimental protocol for determination of fluxes 

We are primarily interested in the high-Rayleigh number limit, where for fixed fluid 
parameters Pr and r the transport properties depend only on the background density 
ratio Rq. This implies the need to use a reasonably tall computational box (L z 3> 1 
in units of d). Our simulations must also permit a large enough number of fingers to 
exist in the horizontal directions to provide good statistical estimates of the turbulent 
fluxes. On the other hand, the domain size should be small enough to suppress any 
secondary large-scale instabilities, which would drive the local density ratio R p away 
from the background Rq and modulate the turbulent fluxes we are trying to measure. 
After a careful study of the outcome of a series of simulations, further detailed below, 
we found that a computational domain of size 5 x 5 x 10 in units of the "fastest growing 
wavelength" (FGW) is a good compromise. The FGW is defined as the wavelength of the 
fastest growing mode of the fing ering ins t ability , and depends on the model parameters 



Pr, r and Rq (see e.g. figure 4 of lSchmittl (|19T9h ) . We find that it is also a good measure 



of the width of the fully nonlinear fingers in the parameter regime considered (see figure 
M and discussion below) . Finally, note that since large-scale perturbations cannot grow 
in these "small-box" simulations, R p = Rq everywhere in the domain. In this section, we 
will equivalently use the notation Rq or R p to denote the density ratio since they are the 
same. 

After selecting a value of the density ratio, calculating the corresponding domain size 
in units of d, and determining an appropriate spectral resolution for the simulations 
(see below, and table [1} , we initialise the calculation with low-amplitude white noise 
perturbations in T and S, and let the system evolve with time. Vertical turbulent fluxes 
of heat and salt, defined as 

F T =< wT >, 

F s =<wS>, (A 2) 

are then computed, where < ... > denotes a volume average over the computational 
domain. Figure [6] shows time series of —Ft and — F$ (note that both fluxes are negative 
for fingering convection), as well as that of their ratio 7, for three values of R p and for a 
fluid with the characteristic properties of salty water (r = 0.01, Pr = 7). After a period 
of exponential growth, the system settles into a statistically stationary finite amplitude 
state in which all plotted quantities fluctuate about well-defined temporal averages, the 
state of homogeneous fingering convection. The turbulent fluxes reported in table [1] are 
measured from temporal means of Ft and F$ once the system is in that state. 



A. 3. The effect of the domain size 

In the last part of this appendix, we estimate the optimal box size for calculating local flux 
laws, while meeting the criteria described above (large enough to provide good statistics, 
small enough to suppress secondary instabilities) . Guidance for our choice can come from 
simulations in the less computationally demanding parameter regime of Pr = 7, r = 1/3. 
Figure [7] shows the averaged Nusselt number and 7 values for simulations from two 
domains at these parameters. The first uses a computational domain of approximately 
5x5x8 fastest-growing wavelengths (FGW) , and the second is much taller, approximately 
6 x 6 x 20 FGW. Comparison of the turbulent fluxes as a function of R p shows strong 
agreement between the two sets of simulations except at the highest and lowest values of 
the density ratio, where some divergence occurs (discussed below). Overall, these results 
suggest that the smallest box-size provides sufficient statistics and a large-enough domain 
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Figure 6. Time series derived from simulations of fingering convection in the heat-salt system 
(Pr = 7, r = 0.01) for the three cases R p — 1.2, R p = 1.5 and R p = 2 (the first and third also 
shown in figure [TJ. a) Fluxes —Ft and —Fs (note that for salt fingering Ft and Fs are both 
negative quantities), b) Turbulent flux ratio 7. 



to yield accurate flux laws at most density ratios. Meanwhile, neither set of simulations 
seem to exhibit any large-scale modes of instability, which also satisfies our requirements. 

In order to understand the difference between the turbulent fluxes measured in the two 
geometries at large density ratio, we examine the geometry of the nonlinear fingers in 
more detail. For an experimental estimate of the height of fingers in a given simulation, 
we look to the autocorrelation of the vertical velocity, averaged over all x and y: 

f(s) = J J J w(x,y : z)w(x : y, z + s) dx dy dz 

We estimate a typical finger height by the distance s beyond which the autocorrelation 
function drops below 0.05, and similarly for the typical finger width. 

Figure [5] shows the results of these calculations for the small (5x5x8 FGW) and 
medium (6 x 6 x 20 FGW) domains, using ten sample points separated by 5000 time 
steps at each value of the background density ratio Rq. By inspection of the results for 
the medium domain, we find that fingers are typically about 0.5-2 FGW in width, and 
4-6 FGW in height, for most values of the density ratio except very close to marginal 
stability. As a result, both box sizes are sufficiently wide to accommodate many fingers. 
The small box is tall enough to contain typically two fingers, while the larger box contains 
approximately five. Both domains seem to be satisfactory except at the largest values of 
Rq (where elevator modes are expected to dominate as Rq — > 1/t), for which the small 
box is too short to contain even a single finger. Large fluctuations in the Nusselt number 
result, shown in figure [9] 

For the purpose of extending mean-field theory to more difficult parameter ranges 
(such as Pr C 1, r < 1 of interest in the astrophysical context), these results provide a 
valuable guide. For most values of i?o, boxes of no more than 5 x 5 x 10 FGW in height 
may be expected to provide robust flux averages, while suppressing large-scale structures 
such as gravity waves that would otherwise swamp the finger field and complicate the 
averaging process. However, as the background density ratio increases toward a com- 
pletely stable stratification, taller boxes are necessary to prevent finger overshoot from 
artificially increasing the measured fluxes. 

One final feature bears comment, namely the sharp increase and sudden decrease 
in finger height variability as Rq increases above 1.6. This transition, apparently not 
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Figure 7. Small and medium box average values of the turbulent flux ratio 7(i?o) and the 
logarithm of the Nusselt number Nu(/?o)- For most _Ro values the calculated small-box averages 
of Nu and 7 closely compare with their medium-box counterparts, even where the finger height 
exceeds the small box height (identifiable in the bottom panel of figure [8|. At large Rq, however, 
the averages diverge, seen in the logarithmic (Nuo — 1) values. 



a function of box size, corresponds to a kink in the j(Rp) curve (see figure [7]), but 
speculation as to its cause is deferred to future work. 
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